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The  Formulation  of  Quantum  Statistical  Mechanics  Based  on  the 
Feynman  Path  Centroid  Density.  III.  Phase  Space  Formalism  and 
Analysis  of  Centroid  Molecular  Dynamics 

Jianshu  Cao  and  Gregory  A.  Voth 

Department  of  Chemistry,  University  of  Pennsylvania,  Philadelphia.  Pennsylvania  19104- 
6323 

ABSTRACT 

The  formuiation  of  quantiim  statistical  mechanics  based  on  the  path  centroid  v-ariable 
in  Fe)rnman  path  integration  is  generalized  to  a  phase  space  perspective,  thereby  including 
the  momentum  as  an  independent  dsmamical  variable.  By  virtue  of  this  approach,  operator 
averages  and  imaginary  time  correlation  functions  can  be  expressed  in  terms  of  an  averaging 
over  the  multidimensional  phase  space  centroid  density.  The  imaginary  time  centroid- 
constrained  correlation  function  matrix  for  the  phase  space  variables  is  then  found  to 
define  the  effective  thermal  width  of  the  phase  sp£u:e  centroid  variable.  These  developments 
;xlso  make  it  possible  to  rigorously  analyze  the  centroid  molecular  dynamics  method  for 
computing  quantum  dynamical  time  correlation  fimctions.  As  a  result,  the  centroid  time 
correlation  function  as  calculated  from  certroid  molecular  dyuEunics  is  shown  to  be  a  well- 
defined  approximation  to  the  exact  Kubo  transformed  position  correlation  function.  This 
-analysis  thereby  clarifies  the  underlying  role  of  the  equilibrium  path  centroid  variable  in  the 
quantum  d>'namical  position  correlation  function  and  provides  a  sound  theoretical  basis 
for  the  centroid  molecular  dynamics  method. 
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1.  Introduction 

In  two  papers*'"  (hereafter  referred  to  as  Paper  I  and  Paper  II)  and  a  Communication.^ 
the  formulation  of  quantum  statistical  mechanics  based  on  the  Feynman  path  centroid 
was  extensively  studied.  These  efforts  originated  from  the  notion  that  the  path  centroid 
variable'*  in  equilibrium  quantum  systems  is  the  most  direct  analogue  to  a  classical  variable 
and  should  therefore  possess  both  formally  interesting  and  computationally  useful  prop¬ 
erties.  Paper  I  developed  the  formal  “cornerstone”  for  the  centroid-based  formulation  of 
equilibrium  properties,  introducing  many  of  the  mathematical  tools  necessary  for  subse¬ 
quent  theoretical  developments.  For  example,  a  formally  exact  theory  was  developed  for 
the  equilibrium  density*  associated  with  the  centroid  variable  (the  so-called  “centroid  den¬ 
sity”).  This  analytical  theory  goes  well  beyond  the  Feynman-Hibbs  variational  theory*  for 
the  partition  function  by  employing  an  infinite-order  diagrammatic  perturbation  expansion 
along  with  resummation  and  renormalization  techniques.  The  antdysis  also  explores  the 
diagrammatic  representation  of  various  approximation  schemes*  ®  for  the  centroid  density 
and  then  systematically  improves  upon  those  schemes.  In  addition  to  the  analytical  the- 
oty  for  the  centroid  density  in  Paper  I,  the  quantum  expressions  for  equilibrium  operator 
averages  and  imaginary  time  correlation  functions  were  reformulated  so  that  the  centroid 
density  occupies  the  role  of  the  underlying  statistical  distribution  function.  Taken  together, 
the  developments  in  Paper  I  represent  a  unified  view  of  equilibrium  quantum  statistical 
mechanics  from  the  centroid  perspective. 

In  Paper  II  and  the  Communication,  the  path  centroid  perspective  was  significantly 
. extended  to  address  the  challenge  of  calculating  quantum  dynamical  time  correlation  func¬ 
tions.  The  most  intriguing  and  promising  result  of  the  dynamical  centroid  analysis  is  a 
method  called  centroid  molecular  dynamics  (centroid  Motivated  by  the  appeal 

of  the  centroid  perspective,  it  was  argued  that  the  quantum  position  correlation  function 
can  be  related  to  a  time  correlation  function  for  the  centroid  variable  with  the  centroid 
■'trajectories  ”  generated  by  classical-like  dynamical  equations  on  an  effective,  temperature- 
dependent.  centroid  potential  energv’  surface.  A  number  of  strategies  were  then  developed 
in  Paper  II  to  compute  the  time  correlation  functions  of  general  coordinate-dependent  oper¬ 
ators.  By  virtue  of  the  centroid  MD  approach,  time  correlation  functions  can.  in  principle. 
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be  computed  for  quantum  many-body  systems  with  a  numerical  effort  that  scales  with 
system  size  in  the  same  manner  as  a  classical  molecular  dynamics  i  MD)  simulation.  In  the 
companion  paper.®  some  numerical  algorithms  are  presented  for  centroid  MD  calculations 
in  general  physical  systems  which  obey  Boltzmann  quantum  statistics. 

The  formal  analysis  in  the  earlier  papers.  was  based  on  the  Feynman  path  integral 
formulation  in  coordinate  space.'  Though  it  can  be  argued  that  momentum  dependent 
quantities  can  in  principle  be  obtained  by  taking  time  derivatives  of  the  position  coordinate, 
complications  arise  due  to  the  time  ordering  of  quantum  operators,  especially  when  mixed 
position  and  momentum  terms  appear  in  those  operators.  A  better  treatment,  therefore, 
requires  a  path  integral  centroid  formulation  in  phase  space  which  would  not  only  generalize 
the  earlier  position  centroid-based  formulation.^"^  but  also  provide  a  remedy  to  the  time 
ordering  problem.  Such  a  formulation  is  contained  in  the  present  paper.  In  a  fashion 
similar  to  Paper  I.  the  imaginary  time-correlated  operators  are  reformulated  as  Gaussian 
averaged  functions  which  are  to  be  averaged  over  the  centroid  phase  space  density.  The  final 
formulas  resemble  those  in  coordinate  space,  although  they  are  defined  here  in  phase  space 
through  a  compact  vector-matrix  notation.  For  completeness,  the  expressions  are  also 
formulated  for  many  degrees-of-fireedom  which  is  simply  an  extension  of  a  one-dimensional 
treatment. 

In  terms  of  the  dynamical  centroid  variable  perspective  for  quantinn  time  correlation 
functions,  a  complete  analysis  of  centroid  MD.  and  its  justification,  again  requires  a  path 
integral  formulation  in  phase  space^  so  that  momentum  can  be  treated  as  an  independent 
dynamical  variable.  Indeed,  the  phase  space  centroid  formulation  in  the  present  paper  has 
its  most  important  application  in  the  analysis  of  centroid  MD  presented  in  Sec.  III.  At  the 
time  when  centroid  MD  was  proposed.^'^  a  relationship  between  the  Fourier  transform  of 
the  real  time  position  correlation  function  and  the  centroid  time  correlation  function  was 
identified.  This  relationship  is  based  on  the  analytical  continuation  of  the  vairiationally 
optimized  local  quadratic  approximation^'^®  to  the  centroid-constrained  imaginary  time 
position  correlation  fimction.^  Since  the  Fourier  relationship  holds  exactly  if  the  centroid 
correlation  function  is  replaced  by  the  Kubo  transformed^  position  correlation  function,  it 
was  speculated  that  the  centroid  correlation  function  must  therefore  be  an  appro.\imation 
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to  the  Kubo  transformed  position  correlation  function.  Tlirough  the  phase  space  centroid 
formulation  developed  herein,  the  relationship  between  the  centroid  correlation  function 
;uid  the  Kubo  transformed  position  correlation  function  is  found  to  be  unique  and  cannot 
be  obtained  unless  the  centroid  is  taken  as  the  dynamical  variable.  In  the  end.  centroid 
MD  turns  out  to  be  exact  for  the  first  two  terms  in  the  Taylor  expansion  of  the  correlation 
function  in  time,  and  the  systematic  error  in  the  method  can  be  identified  at  all  subse> 
quent  orders.  The  deviation  firom  the  exact  quantum  time  correlation  function  is  directly 
proportional  to  the  average  thermal  width  of  the  centroid  “particle” ,  as  well  as  higher- 
order  derivatives  of  the  mean  centroid  force.  Centroid  MD  is  thereby  shown  to  capture 
the  leading  quantum  dynamical  behavior  of  the  position  correlation  function  to  within  a 
well-defined  error. 

As  a  direct  extension  of  the  centroid  MD  theory  in  phase  space,  three  strategies  are 
also  presented  in  the  present  paper  to  calculate  time  correlation  functions  of  general  oper¬ 
ators  which  may  depend  on  both  position  and  momentum.  The  derivations  are  relatively 
straightforward  given  the  similar  results  for  the  correlation  fimctions  of  general  coordi¬ 
nate  dependent  operators  published  in  Paper  II.  Nevertheless,  the  demonstration  that  the 
momentiun  can  be  incorporated  into  the  theory  in  the  seune  fashion  as  the  position  is 
significant  both  from  the  formal  point  of  view  and  for  actual  nmnerical  applications. 

The  present  paper  is  orgainized  as  follows:  In  Sec.  II,  the  equilibriiun  formulatiosn 
of  the  Feynman  path  centroid  density,  operator  averages,  and  imaginary  time  correlation 
functions  in  multidimensional  phase  space  are  described.  In  Sec.  Ill,  the  phase  space 
centroid  perspective  is  then  used  to  analyze  and  more  completely  justify  the  dynamical 
centroid  MD  method  for  the  real  time  position  correlation  fimction.  This  theory  is  then 
extended  in  Sec.  IV  to  formulate  three  approaches  for  computing  general  quantum  time 
correlation  functions.  Some  numerical  examples  are  given  in  Sec.  V.  and  concluding 
remarks  appear  in  Sec.  VI. 

II.  Feynman  Path  Centroid  Formulation  in  Phase  Space 

In  Paper  I.  general  imaginar>’  time  correlation  ftmctions  were  expressed  in  terms  of 
the  coordinate  space  centroid  density  pdqc)  and  the  centroid-constrained  imaginary  time 
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position  correlation  function  CrAT.qd-  Here,  the  formalism  of  Paper  I  will  be  e.xtended 
by  a  phase  space  path  integral  formulation^  so  that  the  momentum  appears  as  an  inde¬ 
pendent  \ariable.  The  development  will  first  involve  the  specification  of  the  phase  space 
centroid  density,  then  the  definition  of  the  centroid-constrained  position,  momentum,  and 
cross-correlation  functions,  and  finally  the  formulation  of  the  expressions  for  equilibrium 
operator  averages  and  general  imaginary  time  correlation  functions.  Unless  specified  oth¬ 
erwise.  the  analysis  will  be  presented  for  an  iV-dimensional  system  with  the  position  and 
momentum  variables  described  by  the  iV-dimensional  column  vectors  q  and  p,  respectively. 
The  generalized  vector  C  is  defined  as  the  collection  of  the  2JV  degrees  of  fireedom  in  phase 
space,  i.e..  (^  =  ip. 

.4.  Phase  Space  Feynman  Path  Centroid  Density 

The  phase  space  centroid  density  can  be  straightforwardly  defined  as 

P.(W  =  /'■••ypar)«(Co-*)eip{-S(C(T)i/B}  , 

where  the  path  centroid  vector  in  phase  space  is  given  by 

1  - 
=  mI  ■ 

The  action  functional  for  the  imaginary  time  phase  space  path  integral  is  given  by^ 

•5[C(7')1  =  j  ■^r)-ip{T)  -^T)^U[C(r)]|  (2.3) 

where  m  is  the  diagonal  particle  mass  matrix  and  ^r)  is  understood  as  the  imaginary 
time  velocity  vector.  The  quantum  partition  function  is  related  to  Eq.  (2.1)  such  that 
Z  =  /  dCcpc(Cc)-  The  phase  space  centroid-constrained  correlation  functions^’^  are  defined 
by  the  2-Y  x  2*V  matrix 


(2.1) 


(2.2) 


'  /•••/  PC(r)5(Cc-^o)exp{-5[C(7)i/h} 

Each  element  of  this  matrix  with  the  indices  “i.y  is  given  explicitly  by  the  2x2  block 
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(2.5) 


f  Cc{A(T)pj(0).gc)  Cr{piiT)qj{0).qc)\ 

\Cc(qi{T)pj{Q).qc)  Cciqi{T)qj{0).qc) ) 

where  ((r)  is  the  quantum  path  fluctuation  with  respect  to  the  centroid  position  Co  i  e.. 

=  Cc  +  C(^)-  The  elements  of  the  centroid-constrained  correlation  ftmction  matrix 
in  Eq.  (2.5)  can  also  be  obtained  by  adding  linear  field  terms  of  the  form  -/(r)  •  |(r) 
and  — ^(r)  •  p(t)  to  the  action  functional  in  Eq.  (2.3)  and  by  then  taking  the  appropriate 
second-order  functional  derivatives  of  -j3Fc[/(T),g(T)]  in  the  limit  /(r)  — *■  0  and  g(r)  — »  0. 
where  —  JFclfir),  ^(r)]  =  In  {pc[/(r),  ^(r)]}  is  the  centroid  free  energy  as  functional  of  the 
two  external  fields.  The  centroid-constrained  correlation  function  matrix  in  Eq.  (2.5)  is 
independent  of  the  momentum  centroid  Pc  if  the  potential  V  is  independent  of  the  momen¬ 
tum  variable.  This  fact  can  be  proven  by  considering  the  Fourier  mode  representation* 
of  the  phase  space  action  functional  in  Eq.  (2.3). 

B.  Operator  Averages 

It  proves  informative  to  first  formulate  the  expression  for  the  equilibrium  average  of 
a  general  operator  in  the  phase  space  centroid  perspective.  This  simple  amalysis  identifies 
the  centroid-constrained  correlation  function  matrix  in  Eqs.  (2.4)  and  (2.5)  as  providing 
the  effective  centroid  “width”  factors  in  phase  space.  In  the  phase  space  path  integral 
perspective.*  the  equilibrium  average  of  an  operator  A  is  given  by  the  expression 


(.4)  =  /•••  /  I>CWA(<(0)1  eip{-SlC(T)j/fi},  (2.6) 

Due  to  the  cyclic  invariance  of  the  trace,  the  operator  can  be  evaluated  at  any  point  along 
the  cyclic  imaginary  time  path  q{T).  The  average  in  Eq.  (2.6)  is  first  re-expressed  so  that 
the  centroid  \’ariable  appears  e.xplicitly  in  the  statistical  averaging,  i.e..^'^° 


(.4)  =  Z-'  j  dCc  y  ■■■  j  PdT)«(C<,-<fo)/l(C(0)!e-5'<<’->l'' 
The  operator  .4(C)  is  then  represented  in  2*'V-dimensionai  Fourier  space,  i.e 


(2.7) 


(2.8) 
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where  .4(  k)  is  the  Fourier  transform  of  the  operator,  k  denotes  the  2iV-dimensional  A:-vector 
{kp,kq).  and  k  ■  s  represents  the  contraction  k  ■  C  =  ^nQi)-  Equation  t2.61 

can  then  be  re-\vritten  as 


—  ^  ^  J ^CcPci^c 


where  the  notation  for  the  action  functional  “S[Cc  +  C(^)]’’  is  understood  to  mean  path 
integration  over  the  fluctuations  around  a  fixed  phase  space  centroid  at  Q.  At  this  point, 
a  cumulant  average  of  the  term  exp[ik  •  C(7’)]  in  the  brackets  in  Eq.  (2.9)  can  be  performed 
in  terms  of  the  path  fluctuation  variable  C(^)-  The  cumulant  average,  for  the  present 
purposes,  is  truncated  at  second-order,  but  it  need  not  be,  and  it  is  assumed  that  the 
variable  ^(r)  exhibits  symmetric  fluctuations  about  the  centroid.  After  performing  the 
inverse  Fourier  transformation  by  integrating  over  k  in  Eq.  (2.7),  the  final  result  for  the 
operator  average  in  the  phase  space  centroid  picture  is  given  by®^ 


A)  ~  (Ac(Cc))p, 


where  the  effective  centroid-dependent  quasiclassical  operator  .4c  is  given  by 


•■lc(Cc)  =  +O)c^(0.<fc) 

= -===i===  /'<fC/l(<c+Oe:tpf-CC,-‘(0.fc)C/2l  .  (2,11) 

y/ det  {27rCr(0,  Qc )l  J  -■ 


Here,  the  vector  \ariabie  C  is  obviously  a  Gaussian  vector  associated  with  a  width  ma¬ 
trix  given  by  Cc(0.  ^c)-  The  above  result  is  a  generalization  of  an  expression  obtained 
previously  for  coordinate-dependent  operators.  This  expression  reveals  the  role  played 
by  the  centroid-constrained  correlation  function  matrix  [Eq.  (2.5)j  in  defining  the  effective 
width  factor  in  phase  space  for  the  centroid  quasiparticle. 
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C.  General  Imaginary  Time  Correlation  Functions 

A  general  imaginary  time  correlation  function  is  defined  as 


Cab^t)  =  (A{t)B(0)) 

=  z-^ J ■■■  J vaT)A[aT)]B[m]exp{-s[ar)]m  •  (2.12) 

where  the  operators  A  and  B  are  general  functions  of  the  2N  variables  (p,  ^  and  the 
imaginary  time  interval  is  bounded  such  that  0  <  r  <  h3.  The  goal  here  is  to  refor¬ 
mulate  the  correlation  function  Cab{t)  in  terms  of  the  phase  space  centroid  density  and 
the  centroid-constrained  correlation  function  matrix  Cc{r.qc)  of  Eq.  (2.4).  Following  the 
analysis  of  Paper  I.  one  first  factorizes  the  expression  to  e.\pose  the  integration  over  the 
centroid  density  and  then  expresses  the  operators  A  and  B  as  Fourier  transforms,  giving 

C«(r)  =  I 

f  •••  f  VQ{t) 

X  - .  ■  )  (2.13) 

/•••/  T>C(r) 

where  k  and  k-C  are  as  defined  following  Eq.  (2.8).  The  centroid-constrained  average  in  the 
l)racketed  term  can  be  expressed  as  a  cumulant  expansion  which,  for  the  present  purpose, 
is  truncated  at  the  second  order.  The  result  is  given  by 


(exp[tA:i  ■  s'(t)  +  zA:2  •  C(0)])c  ^  exp< --  Uri  •  Cc(0.9c)  • 


V 

+  k\  •  Crir,  Qc)  ■  ^2  +  ^2  •  Cc(t.  Qc)  -ki  +  ^2  •  Cc(0,  Qc)  -^2  | ,  (2.14) 


where  the  contractions  can  be  explicitly  written  as 


kn'Cc\r.  Qc)  '  k-m  — 


;  =  i  )  =  1 


Cr(Pi{T)pj{0).qc)  CrXpiir)qj{0).qc) 
Cc{qi{T)Pj{0).qc)  Cr{qi{T)qj{0).qc) 
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.  (2.15) 


In  order  to  carrv  out  the  integrals  over  ki  and  ko  in  Eq.  1 2.131.  a  phase  space  vector 
rotation  is  first  performed  such  that 

Cr(7-)  =  [Cl-r) -t-C(0)]/V2 

C-(r)  =  [C(r)-aO)]/^  .  (2.16) 

and  new  centroid-constrained  correlation  function  matrices  are  defined  as 


=  Cc(0.9c)  -i-  Cc(r.<fc) 

C~(r.qc)  =  Cc(0,?c)  -  Cc(T.qc)  ■  (2.17) 

After  performing  the  integrals  in  Fourier  space,  one  arrives  at  the  following  result  for  the 
imaginary  time  correlation  function  in  the  phase  space  centroid  density  picture: 

(^Ab(t)  -  (PAB('’",G))pe  •  (2.18) 

where  the  centroid-dependent  imaginary  time-correlated  operator  product  PABi^-Q)  is 
defined  as  the  multiple  Gaussian  average 

Pab(t,Q)  =  (.4(d+Ci)5(6  +  C2))  .  •  (2.19a) 

Here,  the  vectors  and  C2  related  to  the  two  Gaussian  vectors  C+  and  C-,  such  that 

0.2  =  i((+±C-)  ■  (2.196) 

v2 

where  c._  and  C-  iiave  Gaussian  width  matrices  given  by  C^(r.  q^)  and  C~(r.  (fc).  respec¬ 
tively.  It  should  be  noted  that  if  the  cumulant  expansion  is  ceirried  out  beyond  quadratic 
order  in  Eq.  (2.14).  a  more  complicated  analytical  form  of  Eq.  r2. 19)  would  be  obtained. 
On  the  positive  side,  however,  the  e.xpression  in  Eq.  (2.19)  simplifies  considerably  if  the 
operators  A  and  B  depend  only  on  the  coordinates  and/or  momenta  of  a  single  particle 
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of  the  system  (i.e..  almost  all  of  the  Gaussian  integrals  can  be  "integrated  out"  of  the 
expression). 

Due  to  the  compact  notation  adopted  in  the  preceding  analysis,  the  final  expressions 
in  Eqs.  (2.18)  and  (2.19)  resemble  those  appearing  in  Paper  I.  Nonetheless,  complications 
arise  when  cross  terms  appear  because  of  the  noncommutation  of  position  and  momentum 
operators.  In  fact,  the  two  off-diagonal  terms  in  the  matrix  Cc{t.  Q)  are  complex  functions 
and,  in  turn,  complex  conjugates  of  each  other.  In  the  expression  for  operator  averages  in 
Eq.  (2.10)  and  for  the  time  correlated  operator  product  in  Eq.  (2.19),  different  operator 
orderings  may  lead  to  different  values.  Though  the  operator  ordering  is  not  explicitly  taken 
into  consideration  in  the  Fourier  transforms  in  Eq.  (2.13).  the  final  expressions  in  terms  of 
the  Gaussian  averages  should  always  be  consistent  with  the  original  choice  of  the  operator 
ordering. 

The  multiple  Gaussiem  average  in  the  centroid-constrained  operator  product  [Eq. 
(2.19)]  may  not  be  easy  to  evaluate,  particularly  if  the  operators  A  ox  B  are  not  poly¬ 
nomials  or  exponentials  in  the  phase  space  variables.  However,  if  one  has  an  analytical 
or  numerical  representation  of  the  centroid  density  Pc(Cc)>  a  quadrature  procedure  can  be 
employed  to  evaluate  Eq.  (2.19)  if  the  Gaussian  averages  cannot  be  evaluated  analytically. 
Alternatively,  if  a  good  approximation  for  the  width  matrix  Ce(0.  gc)  is  in  hand,  a  Monte 
Carlo  procedure  can  be  adopted  for  evaluating  Eq.  (2.18)  simultaneously  with  Eq.  (2.19) 
by  using  a  combined  importance  sampling  function  based  on  both  the  centroid  density  and 
the  multiple  Gaussian  distribution  of  Eq.  (2.19). 

III.  Phase  Space  Analysis  of  Centroid  Molecular  Dynamics 

In  the  earlier  papers.”’^  the  development  of  the  centroid  MD  method  was  guided  in 
part  by  a  variational  analysis  and.  in  part,  by  physical  reasoning  and  intuition.  In  this 
section,  a  more  satisfactory  analysis  and  justification  will  be  provided.  Centroid  MD  is 
essentially  a  classical  MD  method  defined  on  a  quantum  mechanical  effective  potential 
energy  siurface.  WTiile  the  deterministic  nature  of  a  classical-like  MD  algorithm  seems  to 
be  at  odds  with  the  uncertainty  inherent  in  quantum  mechanics,  this  paradox  is  partially 
resolved  in  centroid  MD  by  the  introduction  of  the  centroid  mean  force  which  is  obtained  by 
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first  averaging  over  the  quantum  thermal  path  fluctuations.  Moreover,  the  quasiclassicai 
form  in  centroid  MD  arises  from  a  kind  of  quantum  "pre- averaging"  procedure  which 
is  specifically  suited  for  the  computation  of  the  position  time  correlation  function.  The 
evidence  to  date  strongly  suggests  that  centroid  MD  captures  the  major  featmes  of  the 
ensemble  averaged  correlations  of  the  quantum  dynamical  position  operator.  However, 
two  important  questions  deserve  some  further  attention.  These  are:  ( 1 )  Can  centroid  MD 
be  shown  to  always  give  a  well-defined  approximation  to  quantum  position  correlation 
functions?  And.  (2)  Why  does  the  equilibrium  path  centroid  variable  occupy  such  an 
important  role  in  dynamical  correlation  functions? 

Important  progress  on  the  first  question  can  be  found  in  Paper  II.  In  that  paper,  the 
cinalogy  of  the  equilibrium  path  centroid  variable  to  a  classical  dynamical  degree  of  fi’eedom 
was  motivated  within  the  context  of  a  variational  theory  based  on  effective  quadratic  action 
functionals  in  imaginary  time.^  ®®  Unfortunately,  it  is  not  possible  to  estimate  the  absolute 
accuracy  of  the  centroid  MD  method  using  the  variational  approach  alone.  Fortunately, 
the  hint  of  a  different  approach  is  found  in  the  Appendix  of  Paper  II  where  the  analysis  of 
a  centroid  MD  approximation  for  general  correlation  functions  is  presented  in  the  context 
of  the  Kubo  transformed  time  correlation  ftmction.^  A  similar  approach  will  be  developed 
here  to  demonstrate  that  the  centroid  MD  time  correlation  ftmction  is.  in  fact,  always  a 
well-defined  approximation  to  the  e.xact  Kubo  transformed  position  correlation  function. 
For  notationai  simplicity,  the  analysis  is  restricted  to  a  two-dimensional  phase  space,  but 
it  ctin  be  readily  generalized  through  the  vector-matrix  notation  of  the  previous  section. 

In  centroid  MD.  the  centroid  variable  evolves  according  to  the  classical-like  equation 
of  motion. 


7c(i)  =  Pc{t)fm 
Pc{t)  =  Fc{t)  . 

where  the  instantaneous  centroid  force  is  defined  by“’^ 

Friqc)  =  (/(7c  +  7))c 


(3.1) 


(3.2) 
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Here,  the  /  =  —dVfdq  cind  the  symbol  denotes  the  centroid-constrained  average 
with  the  phase  space  imaginary  time  path  integral  or.  as  an  approximation,  the  eifective 
operator  representation  in  Eq.  (2.11).  With  trajectories  in  hand  obtained  from  the  above 
equations,  the  centroid  MD  time  correlation  function  is  given  by“'^ 


C-(t)  =  (go(t)9c{0))p.  .  (3.3) 

where  the  bracket  with  superscript  Pc  means  that  the  initial  conditions  of  the  centroid 
trajectories  are  averaged  with  the  phase  space  centroid  density  in  Eq.  (2.1).  The  compan¬ 
ion  paper^^  is  devoted  to  the  development  of  several  algorithms  for  evaluating  centroid 
dynamics  [Eq.  (3.1)]  in  general  many- body  systems. 

After  expanding  C*  (t)  as  a  Taylor  series  in  terms  of  t.  and  by  taking  into  consideration 
the  time  reversibility  of  correlation  function,  one  obtains 


where  the  expansion  coefficients  are  expressed  as 


(3.4) 


•  (3.5) 

The  operator  Cc  is  the  classical  Poisson  bracket 


C^A  =  {A.Hc}  (3.6) 

for  a  classical-like  centroid  Hamiltonian  defined  by  He  =  Pc/2  +  Vdqc)  with  the  effective 
centroid  potential  given  by  Vc(qc)  =  -kBTln[pc{qc)/{Tnl2T:h^l3Y^% 

At  this  point,  it  proves  to  be  particularly  informative  to  similarly  analyze  the  Kubo 
transformed  position  correlation  function^ 


+  iT)q{Q))dr 


(3.7) 


which  is  directly  related  to  the  quantum  response  function  (cf.  the  Appendix  of  Paper 


II)  and  the  quantum  dynamical  position,  momentum,  and  cross  correlation  functions. 
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For  example,  by  meJdng  use  of  the  Kubo  transformation  and  the  commutator  relation 
pjm  —  [q,H\lih.  one  can  derive  the  following  Fourier  relations: 

Cqq{uj)  =  {Ti3u/2)[cot)i{hl3uj/2)  +  1]  :  (3.8) 

and 


Cpq  (wir^ )  —  iuJTTlCqq  (uJ ) 

Cpp{u;)  =  u^m"Cqq{uj)  .  (3.9) 

where  p  and  q  stand  for  any  two  elements  of  the  multidimensional  vectors  p  and  q.  Equation 
(3.7)  can  also  be  written  as  a  Taylor  expansion,  giving 


*2n 

m  ■ 

where  the  expansion  coefficients  are  expressed  as 


(3.10) 


1 


(3.11) 


Here.  £  is  a  commutator,  the  quantum  analogue  of  the  Poisson  bracket,  i.e. 


CA  =  .  (3.12) 

After  making  use  of  the  definition  of  the  centroid  variable  in  Eq.  (2.2)  and  the  invariance 
of  trace,  ones  obtains 


=  (gc<£^”?)c)pe  •  (3.13) 

The  centroid  correlation  function  and  the  Kubo  transformed  correlation  function  take 
a  similar  analytical  form  [cf.  Eqs.  (3.5)  ai  ^  (3.13)],  the  difference  being  between  the  terms 
C^qc  and  (£*”<7)c  in  Eqs.  (3.5)  and  (3.13),  respectively.  These  terms  can  be  determined 
(\xpiicitly  and  compared  term-by-term.  The  first  few  terms  are: 
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n  =  0. 


q)c  —  Qc  ■ 


(3.14) 


n  =  1. 


{C^q)c  =  (//m)c  =  Fc/m 

n  =  2. 


(3.15) 


C\)c  =  +  +  +  (3.16) 

for  the  Kubo  correlation  function  (Eq.  (3.13)j,  but  for  the  centroid  correlation  function 
[Eq.  (3.5)], 


£*q,  =  f  (Df;  ,  (3.17, 

m’^  m* 

Therefore,  the  two  leading  terms  of  the  Taylor  expansions  for  the  centroid  and  Kubo  cor¬ 
relation  functions  are  the  same,  the  difference  between  them  beginning  with  the  third  term 
(i.e..  at  order  The  latter  term  will  now  be  taken  as  em  example  to  show  how  to  evaluate 
the  leading  correction  term  to  the  centroid  correlation  function  »i.e..  to  demonstrate  that 
the  centroid  correlation  fimction  is  a  well-defined  approximation  to  the  Kubo  correlation 
function).  The  Gaussian  representation  of  operators  in  phase  space  [Eq.  (2.11)]  proves  to 
be  useful,  though  not  essential,  in  the  following  analysis. 

It  is  first  noted  that  the  centroid  average  of  a  product  of  operators  can  be  written  as 
a  product  of  the  centroid  averages  of  operators  and  a  leading  correction  term.  i.e.. 

{AB)c  ^  ArBc  +  AlB'^Cr.  .  (3.18) 

where  terms  higher  order  in  the  phase  space  width  factor  Cc  [Eq.  \2.4)]  have  been  omitted. 
By  virtue  of  the  cyclic  property  of  the  trace,  it  can  be  shown  in  general  that 

C^(p(0)9(0).(?J+a{9{0)p(0),gJ  =  0  .  (3.19) 
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It  is  also  noted  that  the  application  of  the  commutator  £  of  Eq.  (3.12)  four  rimes  in 
Eq.  (3.16)  leads  to  a  symmetrized  arrangement  of  momentum  and  coordinate  operators. 
Combining  this  fact  with  Eqs.  (3.18)  and  (3.19).  it  is  seen  that  there  will  exist  no  centroid- 
(  onstrained  correlation  functions  mixing  momentmn  and  coordinate  operators— at  least  to 
the  order  of  the  leading  correction  term.  Thus.  Eq.  (3.16)  becomes 

=  +  .  (3.20) 

where  stands  for  the  centroid-constrained  operator  average  =  (/^”^)c.  Equation 
(3.18)  is  useful  when  comparing  the  second  term  on  the  right-hand-side  of  Eq.  (3.20)  with 
the  second  term  of  Eq.  (3.17).  i.e..  ./'^V)c  —  The  correction  factor  for  this  term 

contains  the  width  factor  and  second-  or  higher-order  derivatives  of  the  force.  The  focus 
of  the  analysis  can  now  shift  to  a  comparison  of  the  first  term  on  the  right-hand-sides  of 
Eqs.  (3.17)  and  (3.20). 

Because  the  centroid  force  Fc  is  a  function  and  /  is  a  quantum  operator,  the  n-th 
order  derivative  of  the  centroid  force.  and  the  centroid-constrained  average  of  the 
n-th  order  deri^ative  of  the  force,  ,  are  different  even  though  the  centroid  force  is  equal 
to  the  centroid-constrained  average  of  the  force,  i.e.,  Fc  =  (/)c-  .Application  of  the  chain 
rule  to  reveals  the  difference  between  the  force  derivative  terms  as 

=  -5/i‘'""’Cj‘’(«(0)9(0).?,)  +  ---  .  (3.21) 

—  where  the  deri\ative  of  the  width  factor  appears  here  instead  of  the  width  factor  itself.  This 
observation  is  particularly  significant  because  it  shows  that  an  e.xpression  in  terms  of  the 
centroid  force  or  its  derivatives  agrees  with  its  quemtiun  counterpart  to  all  orders  in 
the  width  factor  Gc,  with  corrections  coming  only  in  the  spatial  derivatives  of  that  factor. 
From  Eqs.  (3.5)  and  (3.13),  it  is  seen  that  the  latter  correction  is  then  to  be  averaged 
over  the  centroid  density,  so  large  deviations  will  occur  only  if  the  width  experiences  large 
liuctuations  which  persist  in  the  average  sense. 

Returning  to  the  first  term  in  the  right-hand-side  of  Eq.  (3.20),  the  quantmn  fluc¬ 
tuations  in  momentum  will  contribute  a  further  deviation  from  the  similar  term  in  the 
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centroid  correlation  function  [Eq.  (3.17)].  The  difference  between  the  two  terms  for  all 
powers  of  n  is  given  by 


^Pc  =  (pVc-Pc  =  ^Pc  '^Ccip{0)p{0}.qc)  +  ---  .  (3.22) 

where  Cc(p(0)p{0),  Qc)  is  the  width  factor  for  the  momentum  path  fluctuations.  The  terms 
associated  with  this  correction  have  a  value  of  n  no  smaller  than  2.  The  average  of  the 
term  phase  space  centroid  density  will  simply  factorize  and  give  a  constant 

since  the  distribution  of  the  momentum  centroid  is  the  classical  Boltzmann  distribution 
imd  independent  of  the  higher-order  path  Fourier  modes.^*^ 

Taking  into  account  ail  of  the  preceding  considerations  and  generalizing  them  to  terms 
of  higher  order,  one  can  sum  up  in  general  the  difference  between  the  Kubo  transformed 
position  correlation  hmction  [Eq.  (3.11)|  and  the  centroid  MD  correlation  function  [Eq. 
(3.3)  to  give 


^  *2n 

=  ,  (3.23) 

with  A  being  the  difference  between  the  two  sets  of  coefficients,  given  by 

^12n|  ^  ^I2ni  _  ^  •  (3.24) 

I 

where  at  least  one  n;  is  no  smaller  than  2  and  terms  which  involve  the  spatial  derivatives  of 
the  width  factor  Cc  have  been  neglected.  The  above  result  can  be  confirmed  by  dimensional 
..-—analysis. 

The  preceding  analysis  confirms  that  the  apparent  success  of  centroid  MD  is  by  no 
means  incidental  and  is.  in  fact,  both  physically  and  mathematically  understandable.  A 
few  comments  on  the  implications  and  significance  of  the  analysis  are  as  follows: 

(a)  At  first  glance.  Eq.  (3.24)  confirms  two  properties  of  centroid  MD  which  are  already 
obvious."-^  Namely,  the  method  is  clearly  exact  if  the  potential  contains  no  global  anhar- 
monicity  (i.e..  =  0)  or  if  the  system  is  near  the  classical  limit  (i.e..  Cc  is  small). 

However,  upon  closer  inspection  Eq.  (3.24)  also  reveals  why  centroid  MD  is  a  good  ap¬ 
proximation  in  other,  less  straightforward,  situations.  To  be  specific,  the  centroid  force  Fc 
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[Eq.  (3.2)]  at  the  centroid  position  Qc  is  computed  by  averaging  the  classical  force  over  the 
imaginary  time  Feynman  paths  around  the  centroid.  This  averaging  tends  to  smooth  out 
rather  dramatically  any  "kinks"  in  the  potential  energy  function,  at  least  over  the  length- 
scale  of  the  particle  s  thermal  width.  The  more  "quantum"  the  particle,  the  more  the 
averaging  occurs.  This  behavior  is  very  important  and  is  one  reason  an  effective  quadratic 
action  functional  can  be  so  accurate  in  describing,  e.g.,  the  equilibrium  properties  of  the 
polaron^^  and  the  hydrated  electron. Given  that  such  local  "smoothing”  occurs,  the 
higher-order  derivatives  of  the  centroid  force  in  Eq.  (3.24)  will  tend  to  be  small,  giving  a 
small  correction  factor  in  Eq.  (3.24)  even  though  the  quantmn  thermal  width  Cc  may  be 
sizable.  In  order  for  centroid  MD  to  become  inaccurate,  the  preceding  argument  suggests 
that  the  locally  nonlinear  features  of  the  effective  centroid  potentiad  (i.e..  the  effective  an- 
harmonic  terms)  must  be  large.  Such  behavior  occurs  if  the  local  anharmonicities  in  the 
physical  potential  near  the  centroid  are  large  in  the  average  sense  and  persist  on  a  length- 
scale  greater  than  the  thermal  width  factor  of  the  centroid  quasiparticle.  Moreover,  since 
any  correction  terms  are  then  averaged  in  Eq.  (3.24)  over  the  centroid  distribution,  the 
effective  anharmonic  behavior  must  also  be  present  to  a  significant  degree  in  the  regions 
of  greatest  centroid  density.  Interestingly,  as  the  system  becomes  more  classical  the  effec¬ 
tive  anharmonicities  (i.e.,  the  higher  derivatives  of  the  centroid  force)  will  indeed  become 
larger,  but  this  effect  will  be  compensated  for  by  a  greater  reduction  in  the  thermal  width 
factor  Cc  in  Eq.  (3.24).  (Recall  that  centroid  MD  always  yields  the  exeict  classical  limit.) 

(b)  The  definition  of  the  centroid  force  captures  the  major  contribution  of  quantum  fluc¬ 
tuations  and  predicts  the  right  quantum  dynamics  to  within  a  tolerance  proportional  to 
the  averaged  higher-order  derivatives  of  the  centroid  force  and  the  centroid  thermal  width 
fawtor  Cc.  Though  most  semiclassical  approximations  are  expanded  in  terms  of  h.  the 
centroid  MD  correlation  function  already  contains  terms  which  are  infinite  order  in  h.  The 
leading  corrections  to  the  centroid  dynamics  depend  on  the  thermal  width  factor  for  the 
potential  at  hand  (i.e.,  not  just  the  free  particle  width  which  is  of  order  h^).  In  the  nearly 
classical  limit,  the  centroid  dynamics  approximates  the  quanttun  dynamics  to  the  order  of 

instead  of  h  as  in  the  WKB  approximation. 

(c)  In  general,  the  Kubo  transformed  position  correlation  fimction  [Eq.  (3.7)]  is  a  weil- 
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defined  quantum  quantity  which  is  an  ideal  candidate  for  any  t>'pe  of  a  classical-like  ap¬ 
proximation.  This  property  arises  because  the  integration  over  the  imaginary  time  r 
eliminates  the  imaginary  part  of  the  correlation  function  and  also  averages  out  certain 
quantum  fluctuations. 

(d)  To  shed  some  light  on  the  unique  role  of  the  centroid  in  formulating  a  classical-like 
approach  to  the  position  correlation  function,  one  can  apply  the  Taylor  expansion  to  the 
symmetrized  position  correlation  function  C(t),  giving 


<?W  =  i([«(t)-«(0)W  =  ,  (3.25) 

“  n=tt) 

where  I*  •  •]+  is  the  anticommutator.  If  oue  attempts  to  carry  out  a  term-by-term  analysis 
as  for  the  Kubo  transformed  position  correlation  function,  difficulties  immediately  arise 
with  the  first  few  terms  because  it  is  not  clear  how  to  define  the  momentum  distribution 
in  the  general  case.^^  (Note  again  tnat  the  centroid  momentum  distribution  is  simply  the 
Boltzmann  distribution.^*^)  Consequently,  it  is  not  clear  from  this  analysis  of  Eq.  (3.25) 
how  to  define  an  accurate  classical  MD-like  algorithm  to  compute  the  symmetrized  quan¬ 
tum  correlation  function  for  general  potentisus.  Focusing  instead  on  the  Kubo  transformed 
position  correlation  function  [Eq.  (3.7)]  reveals  the  factorization  of  the  centroid  \*ariable 
which,  in  turn,  leads  one  to  the  factorization  of  the  centroid  constrmned  average  in  Eq. 
(3.13).  The  subsequent  identification  of  the  centroid  force  in  the  Taylor  expansion  terms 
w>.-of  the  correlation  function  supports  the  conclusion  that  the  centroid  variable  can  indeed 
be  viewed  as  a  dynamical  variable  at  a  well-defined  level  of  approximation. 

(e)  To  impro^•e  its  accuracy,  it  seems  certain  that  centroid  MD  should  be  augmented  by 
im  additional  quantum  factor.  Because  the  correction  to  the  centroid  force  begins  at  the 
H  term  in  the  Taylor  series  expansion,  such  a  term  will  not  add  linearly  to  the  determin¬ 
istic  centroid  force,  but  it  might  instead  be  constructed  as  some  kind  of  time  convolution 
reflecting  the  non-deterministic  nature  of  quantum  mechanics.  .Apparently,  this  "quantum 
menKHT  function*  would  depend  locally  on  the  width  factor  and  the  anharmonicity  and 
stil^  d  a  time-reversible  dynamics.  Of  course,  this  argument  is  ptirely  speculation. 
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rV.  General  Time  Correlation  Functions 

In  Paper  II.  three  different  centroid-based  approaches  were  proposed  for  calculating 
time  correlation  functions  consisting  of  operators  which  depend  on  the  coordinate  \'ari- 
jibles.  In  this  section,  these  strategies  are  modified  to  allow  the  computation  of  correlation 
functions  of  general  operators  which  may  depend  on  both  position  and  momentum.^'^  The 
phase  space  centroid  perspective,  developed  in  Sec.  I.  will  be  employed  to  extend  the 
formulations  of  Paper  II  to  general  operators.  The  results  in  the  present  case  will  also  be 
given  for  a  multidimensional  phase  space  using  the  notation  in  accord  with  that  introduced 
in  Sec.  11.  In  addition,  the  companion  paper®  describes  several  algorithms  for  computing 
centroid  MD  trajectories. 

A.  Analytical  Continuation  of  Centroid- Constrained  Correlation  Functions 

One  of  primary  results  from  Sec.  II  is  the  expression  for  general  imaginary  time 
phase  space  correlation  functions  in  the  centroid-based  perspective  [Eqs.  (2.12)-(2.19)]. 
In  principle,  the  double  Gaussian  average  in  Eq.  (2.19)  can  be  performed  for  any  func¬ 
tional  form  and  the  resulting  expression  will  then  involve  an  average  of  functions  which 
depend  on  Cc(r,  qc)  over  the  normalized  centroid  density.  At  that  point,  the  centroid- 
constrained  propagator  can  be  analytically  continued  into  real  time  (r  — ►  it)  to  yield  an 
approximation  to  the  rezd  time  correlation  function  Within  the  framework  of  the 

approximate  locally  optimized  effective  harmonic  theory.^  ®  the  real  time  version  of  the 
centroid-constrained  correlation  function  matrix  can  be  obtained  by  replacing  r  with  it 
,,..and  the  resulting  expressions  used  in  Eqs,  (2.18)  and  (2.19). 

.\s  was  pointed  out  in  Paper  II.  this  tmalytical  continuation  method  may  not  be 
completely  satisfactory  for  two  reasons:  (1)  The  effective  harmonic  version  of  Eq.  (2.4) 
is  dynamically  accurate  only  for  relatively  short  times.  The  anharmonicities  in  the  real 
potential  wiU  cause  the  analytically  continued  effective  quadratic  correlation  function  to 
deviate  from  the  exact  behavior  at  long  times  even  in  the  classical  limit:  (2)  The  correlated 
operator  representation  in  Eq.  (2.18)  is  expressed  at  the  level  of  a  second-order  cumulant 
expansion.  Though  tliis  approximation  may  be  an  excellent  one  for  imaginary  time  cal¬ 
culations.  approximate  re£d  time  correlation  functions  can  be  more  sensitive  to  nonlinear 
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interactions  and  thus  less  stable  in  their  behavior  at  long  times. 


B.  Camuiant  Expansion  Combined  with  Centroid  AID 

As  has  been  shown  in  Sec.  III.  the  centroid  MD  correlation  function  lEq.  (3.3)]  gener¬ 
ally  provides  an  accurate  representation  of  the  exact  Kubo  transformed  position  correlation 
function  which,  in  turn,  yields  the  real  time  correlation  functions  (C(f)C(O))  through  the 
Fourier  relations  in  Eqs.  (3.8)  and  (3.9)  after  replacing  with  C*(lo).  Another  strategy 
for  computing  general  correlation  functions  is  therefore  to  first  introduce  the  phase  space 
variable  correlation  functions  directly  into  the  expression  for  the  general  correlation  func¬ 
tion  {A{t)B{0))  and  to  then  use  centroid  MD  to  calculate  (C(t)C(O))  in  that  expression. 
Though  many  of  the  following  expressions  in  this  approach  are  similar  to  those  given  in 
Paper  II.  the  derivation  in  phase  space  is  included  here  for  completeness. 

To  begin  the  derivation,  one  considers  the  general  imaginary  time  correlation  function 
Cab(t)  =  (A(r)B(O))  and  expresses  it  as 

CAB{r)  =  j  J  {exp[iki '  ^{t)  + 'k2  •  CiO)])  ,  (4.1) 

where  the  canonical  path  integrtd  average  in  phase  space  is  explicitly  given  by 

(...)  =  /•••/PC(r)  (•••)  exp{-5fC(r)]/^} 

/  •  •  •  /  '^Mexp{-S[C{T)]/h} 

The  cumulant  average  of  the  exponentitd  term  in  Eq.  (4.1)  can  be  performed  and.  for 
.-Simplicity,  truncated  at  second  order,  giving 


{exp[iki  ■  ;{t)  -I-  ik2  •  C(0)1) 


expl  iki  ■  (C)  +  ik2  •  (0  -  - 


ki  ’  Cf,{0)  •  ki  -r  ki  •  Cs(t)  •  ko 


^2  ’  '  ki  ~  k2  •  C^(0)  •  ko 


.  (4.3) 


where  the  imaginary  time  fluctuation  correlation  functions  constitute  a  Hermitian  matrix, 
defined  by  the  elements 
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[C6(r)]ij 


{(*=Pi(T)<!*>j(0))  {6pi{T)dQj{0)) 

{dqi{T)6pj{Q))  {6qiir)6qj{0)) 


(4.4) 


where  6C  =  C  -  C)-  In  order  to  perform  the  integrals  over  A:i  and  ko  in  Eq.  (,4.1).  new 
imaginary  time  correlation  function  matrices  are  defined  as 

C±(r)  =  C«(0)±C«(r)  .  (4.5) 

After  performing  the  A:<integrals  in  Eq.  (4.3),  the  expression  for  the  general  imaginary  time 
correlation  function  is  given  by  the  double-Gaussian  average 

C.,b(t)  ;=  (.■l«d+Ci)B((d+C2))^^^.  .  (4.6) 

5  *  5 

where  the  vectors  4'i  and  C2  are  related  to  the  two  Gaussiem  vectors  C+  and  C-  such  that 

d.2  =  ■i(d+±C-)  .  (4.7) 

The  Gatissian  vectors  C+  and  C-  have  width  matrices  given  by  CJ{t)  and  C7(t),  respec¬ 
tively,  as  defined  in  Eq.  (4.5). 

The  imaginary  time  expressions  for  C*(r)  can  now  be  analytically  continued  via  the 
inverse  Wick  rotation  r  —  it.  The  resulting  approximate  expression  for  the  real  time 
correlation  function  Cab  it)  is  given  by 

CAB(t)  ^  (A((0+Ci)S((C)+C2))^^,,^_,  .  (4.8) 

\  /C7(t),c,  (t) 

where  the  real  time-dependent  Gaussijin  width  matrices  are  given  by 

Cf[t)  =  Q(0)±Q(t)  .  (4.9) 

The  correlation  function  elements  of  Eq.  (4.9)  can  be  calculated  using  the  centroid  MD 
position  correlation  function  C*(f)  defined  in  Eq.  (3.3)  and  the  Fourier  transform  relations 
in  Eqs.  (3.8)  and  (3.9)  with  miu)  replaced  by  C'*(u;).  Equation  (4.8),  with  Eq.  (4.9).  is  the 
central  result  of  this  subsection.  It  should  be  noted  that  Eq.  (4.8)  simplifies  considerably 
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if  the  operators  .-1  and  B  depend  on  only  one  phase  space  coordinate  (i.e..  most  of  the 
Gaussian  integrals  can  be  integrated  out  of  the  expression). 

C.  Centroid  MD  with  Semiciassicai  Operators 

A  third  algorithm  is  discussed  in  this  subsection  for  calculating  general  time  correlation 
ftmctions  with  centroid  MD.  This  procedure,  which  is  more  approximate  than  the  previous 
ones,  assumes  a  semiciassicai  representation  of  the  quanttim  operators,  the  goal  here  being 
one  of  computational  utility. 

The  method  called  “centroid  MD  with  semiciassicai  operators" ,  introduced  in  Paper 
II.  centers  around  the  computation  of  the  real  time  correlation  fimction  C*^ff{t),  given  by 

=  (Ae(f)5e(0)>p,  .  (4.10) 

where  the  initial  condition  averaging  is  performed  with  the  (normalized)  phase  space  cen¬ 
troid  density  defined  in  Eq.  (2.1).  The  semiciassicai  operators  Ociqdt))  in  Eq.  (4.10)  are 
given  by  the  time-dependent  analog  of  Eq.  (2.11),  i.e.. 

<^c(^(<))  =  (0(G(^)  +C))Cc(0,fe(t))  (4.11) 

Here.  Cc(0  is  the  phase  space  centroid  trajectory  which  obeys  the  centroid  MD  equation 
of  motion  in  Eq.  (3.3),  and  the  time  dependent  Gaussian  width  matrix  0^(0,  qdt))  for  the 
vector  C  is  given  by  the  centroid-constrained  correlation  function  matrix  in  Eq.  (2.4)  with 
the  position  centroid  qo  located  at  qdt). 

As  shown  in  the  Appendix  of  Paper  II,  the  general  centroid  correlation  ftmction  in 
Eq.  (4.10)  is  an  approximation  to  the  Kubo  transformed  version  of  the  exact  correlation 
function  CAB(t).  Therefore,  in  order  to  calculate  CAB(t)  one  makes  use  of  the  Fourier 
relationship" 


OabH  ^  (n3u/2)lcotli(n0u;/2)  +  1]  (4.12) 

The  expression  in  Eq.  (4.10)  is  intended  to  maximize  the  utility  of  the  centroid  MD  method 
in  a  transparent,  though  approximate,  fashion  for  general  correlation  functions.  The  reader 
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is  referred  to  Paper  II  for  further  analysis  and  comments. 


V.  Numerical  Examples 

In  this  section,  the  accuracy  of  the  centroid  MD  formulation  is  tested  for  correlation 
functions  other  than  the  position  correlation  function.  In  these  studies,  a  nonlinear  oscil¬ 
lator  model  is  employed  which  is  the  same  as  one  studied  in  Paper  II.  The  potential  for 
this  model  is  given  by 

ViQ)  =  ^q'+cq^+9Q^  (5.1) 

with  the  parameters  c  =  0.10.  g  =  0.01.  h  =  1.0.  and  m  =  1.0.  In  these  units,  the  inverse 
temperature  d  is  given  as  values  of  the  dimensionless  parameter  dhu.  The  single  minimum 
of  the  potential  in  Eq.  (5.1)  is  located  at  g  =  0.  The  cubic  anharmonicity  is  operational 
for  small  deviations  from  the  minimum,  while  the  quartic  anharmonicity  influences  the 
larger  deviations  from  the  minimum.  At  low  temperatmres.  the  cubic  anharmonicity  is  the 
dominant  perturbation.  The  energy  gap  between  the  ground  and  first  excited  vibrational 
state  for  the  potential  in  Eq.  (5.1)  is  shifted  to  the  red  by  6%  compared  with  the  harmonic 
limit  of  the  energy  spectrum.  Such  an  anharmonic  shift  is  equivalent  to  a  (rather  large) 
shift  of  180  cm'^  for  a  3000  cm'^  C-H  stretching  mode  or  a  60  cm~^  shift  of  a  1000  cm~^ 
C-C  stretching  mode.  A  temperature  of  d  =  10  was  employed  in  the  calculations  which 
is  equivalent  to  a  C=C  double  bond  at  300  K. 

The  exact  correlation  functions  for  the  potential  in  Eq.  (5.1)  were  obtained  by  di- 
— agonalizing  the  nonlinear  Hamiltonian  in  a  harmonic  oscillator  basis  and  employing  one 
hundred  of  the  eigenstates  in  the  dynamics  calculations.  The  centroid  forces  and  centroid 
potential  in  the  centroid  MD  calculations  were  calculated  from  the  optimized  harmonic 
reference  centroid  density.^-®  This  approximation  was  shown  in  Paper  I  to  be  an  extremely 
good  representation  of  the  exact  result.  The  centroid  potential  and  forces  were  interpo¬ 
lated  from  a  1000-point  grid  within  the  range  (—10. 10] .  The  centroid  MD  initial  positions 
were  generated  by  Metropolis  importance  sampling  from  the  position  centroid  distribution, 
while  the  initial  momenta  were  directly  sampled  from  the  Gaussian  centroid  momentum 
distribution  function.  The  evolution  of  10'’  centroid  trajectories  were  then  calculated  us- 
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ing  the  leap-frog  algorithm  with  a  time  step  of  0.05.  The  classical  MD  simulations  were 
performed  in  the  same  fashion  except  that  the  real  potential  and  force  were  used  instead 
of  centroid  quantities. 

.4.  Momentum  Correlation  Functions 

In  Fig.  1.  the  momentum  correlation  function  is  shown  for  the  nonlinear  oscillator 
in  Eq.  (5.1)  at  J  =  10.  The  solid  circles  are  the  exact  quantum  results,  while  the  solid 
line  is  the  centroid  MD  result.  The  latter  method  is  clearly  very  accurate.  It  should  be 
noted  that  results  of  similar  accuracy  were  reported  in  Paper  II  for  the  position  correlation 
function. 

B.  General  Correlation  Functions 

In  order  to  test  the  methods  outlined  in  Sec.  IV  for  calculating  general  correlation 
ftmctions  in  the  phase  space  centroid  perspective,  the  correlation  function  (p^(t)p^(O))  was 
calculated  for  the  nonlinear  potential  defined  in  Eq.  (5.1).  This  correlation  fimction  is  a 
serious  test  for  the  approximate  methods  because  of  its  nonlinearity  and  the  fact  that  its 
classical  amplitude  is  almost  completely  negligible  at  lower  temperatures.  In  Fig.  2.  the 
results  for  a  temperature  of  /3  =  10  are  shown.  The  soUd  circles  are  the  exact  quantum 
results,  rhe  solid  line  is  the  cumulant  expansion  with  centroid  MD  theory  of  Sec.  IVB. 
the  dashed  line  is  the  centroid  MD  with  semiclassicai  operators  result  from  Sec.  IVC, 
and  the  dot-dashed  line  is  the  classical  MD  result.  The  correlation  function  from  the 
_^xumuiant  expansion  theory  probably  obtains  the  best  agreement  with  the  exact  result, 
although  there  is  too  much  structure  and  synunetry  in  the  oscillations.  The  centroid  MD 
with  semiclassicai  operators  method  is  also  accurate,  but  it  does  not  seem  to  reproduce 
the  higher  frequency  oscillation.  By  comparison,  the  classical  MD  result  is  extremely 
inaccurate.  Qualitatively  similar  results  were  obtained  in  Paper  II  for  the  correlation 
function  (q^{t)q-"{Q)).  A  discussion  of  the  strengths  and  weaknesses  of  the  cumulant  with 
centroid  MD  method  was  provided  in  Paper  II.  and  this  discussion  is  also  applicable  to 
the  phase  space  centroid  formulation  of  this  method. 

In  Fig.  3.  the  correlation  function  (A(f)B(O)).  where  .4  =  pq  and  B  =  qp,  is  shown  for 
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the  nonlinear  potential  in  Eq.  (5. 1 )  at  d  =  10.  The  symbols  and  lines  are  the  same  as  in  Fig. 
2.  This  correlation  function  presents  another  serious  test  of  the  various  methods  because  of 
the  non-commutation  of  the  position  and  momentum  operators.  Again,  the  classical  MD 
result  is  extremely  inaccurate  for  this  low  temperature  correlation  fimction.  The  centroid 
MD  with  semiclassictii  operators  method  does  not  reproduce  the  amplitude  and  negative 
values  of  this  correlation  function.  This  feature  of  the  latter  method  arises  because  the 
correlation  of  the  two  operators  at  different  times  is  ignored  when  the  Gaussian  averages 
are  performed.  Consequently,  the  semiclassicai  operator  approximation  underestimates  the 
real  time  interference  of  the  two  quantum  operators.  On  the  other  hand,  the  accuracy  of 
the  centroid  MD  with  semiclassicai  operator  method  is  superior  to  the  classical  calculation. 
Apparently,  only  the  cumuiant  method  can  describe  the  quantum  interference  effects  for 
this  correlation  function,  and  it  appears  to  do  so  quite  well. 

Though  the  centroid  MD  results  are  quite  encouraging  and  far  superior  to  classical 
MD.  it  seems  evident  that  none  of  the  centroid  MD  approaches  developed  in  Sec.  IV  for 
general  correlation  functions  are  completely  satisfactory  under  all  circumstances.  Future 
research  will  be  devoted  to  this  important  issue. 

VI.  Concluding  Remarks 

In  the  present  paper,  the  formulation  of  quantum  statistical  mechanics  b.  sed  on  the 
path  centroid  \'ariable  in  Feynman  path  integration  has  been  generalized  to  a  phase  space 
perspective.  By  virtue  of  this  perspective,  one  can  express  operator  averages  and  imaginary 
—  time  correlation  functions  in  terms  of  a  classical-like  averaging  over  the  miiitidimensional 
phase  space  centroid  density.  .4n  imaginary  time  centroid-constrained  correlation  function 
matrix  for  the  phase  space  variables  is  seen  to  provide  the  effective  width  factors  for  the 
phase  space  centroid  variables.  The  most  significant  aspect  of  the  phase  space  analysis  is 
that  it  facilitates  a  rigorous  analysis  and  justification  of  the  centroid  molecular  dynam¬ 
ics  method  for  computing  quantum  dynamical  time  correlation  functions."’^  Specifically, 
the  centroid  time  correlation  function  calculated  with  centroid  MD  is  shown  to  be  a  well- 
defined  approximation  to  the  exact  Kubo  transformed  position  correlation  function.  This 
analysis  therebv  reveals  the  important  and  not  completely  obvious  connection  between 
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the  equilibrium  path  centroid  variable  and  the  quantum  dynamical  position  correlation 
function.  Several  strategies  were  then  be  developed  for  using  centroid  MD  in  the  compu¬ 
tation  of  general  time  correlation  functions  of  quEintum  operators  which  depend  on  both 
position  and  momentum.  In  addition,  the  companion  paper®  describes  several  algorithms 
for  computing  centroid  dynamics  in  realistic  many-body  systems.  Future  publications  will 
be  devoted  to  the  application  of  centroid  MD  in  the  simulation  of  rezd  systems,  as  well  as 
to  the  continuation  of  its  formal  development. 
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Figure  Captions 


Fig.  1:  A  plot  of  the  real  time  momentum  autocorrelation  function  for  the  nonlinear  oscillator 
described  in  Eq.  (5.1)  at  a  temperature  of  =  10.  The  solid  circles  are  the  exact 
quantxun  results,  while  the  solid  line  is  the  centroid  MD  result.^"* 

Fig.  2:  A  plot  of  the  correlation  function  {p^(t)p^(0))  for  the  nonlinear  potential  in  Eq.  (5.1) 
at  a  temperature  of  3  =  10.  The  solid  circles  are  the  exact  quantum  results,  the  solid 
line  is  the  cumulant  expansion  with  centroid  MD  theory  of  Sec.  IVB,  the  dashed  line  is 
the  centroid  MD  with  semiclassical  operators  method  of  Sec.  IVC.  and  the  dot-dashed 
line  is  the  classical  MD  result. 

Fig.  3:  A  plot  of  the  correlation  function  {A{t)B{0)),  where  A  =  pq  and  B  =  qp,  for  the 
nonlinear  potential  in  Eq.  (5.1)  at  a  temperature  of  13  =  10.  The  solid  circles  are 
the  exact  quantum  results,  the  solid  line  is  the  cumulant  expansion  with  centroid  MD 
theory  of  Sec.  IVB,  the  dashed  line  is  the  centroid  MD  with  semiclassical  operators 
result  from  Sec.  IVC,  and  the  dot-dashed  line  is  the  classical  MD  result. 
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